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Abstract 

A three-dimensional extension of the structural default model with 
firms' values driven by correlated diffusion processes is presented. Green's 
function based semi-analytical methods for solving the forward calibration 
problem and backward pricing problem are developed. These methods are 
used to analyze bilateral counterparty risk for credit default swaps and 
evaluate the corresponding credit and debt value adjustments. It is shown 
that in many realistic cases these value adjustments can be surprisingly 
large. 



1 Introduction 



1.1 Motivation 



The recent turmoil in financial markets has profoundly changed their modus 
operandi. Credit trading in general, and correlation trading in particular, un- 
derwent important transformations. Standardization of credit default swaps 
(CDSs) and the development of clearing houses for their trading are just two 
examples of recent changes aimed at a more transparent setup in the credit 
market. At the same time, trading volumes for bespoke tranches of collateral- 
ized debt obligations (CDOs) have shrunk significantly compared to the peak in 
2007; while more complex structures such as CDOs-Squared have almost disap- 
peared. The focus has shifted from more complicated products, towards simpler 
products, such as credit indices, collateralized CDSs, funded single name credit- 
linked notes (CLNs), CDSs collateralized by risky bonds and other products, 
for which risks are somewhat easier to understand, model, and mitigate. More 
details can be found in several recent books, including Berd 
et al.| |201l| , |Gregory| [2011| , |Lipton and Rennie| [2011] 



2010 , Bielecki 



As a result of the financial crisis, the need for proper accounting of counter- 
party risk in the valuation of over-the-counter (OTC) derivatives has become 
paramount. This has happened due to the fact that some protection sellers, 



such as mono-line insurers and investment banks, have experienced sharply ele- 
vated default probabilities or even default events, the case of Lehman Brothers 
being the prime example. Counterparty credit risk can be defined as the risk 
of a party to a financial contract defaulting prior to the contract's expiration 
and not fulfilling all of its obligations. This risk can be mitigated by collat- 
eralizing the corresponding contract or moving it to an exchange. However, 
in some cases this is not possible, and many OTC contracts are privately ne- 
gotiated between counterparties and subject to counterparty risk. Since both 
parties to a particular contract can default, one needs to account for both credit 
and debt value adjustments. The valuation of OTC products poses a common 
problem: companies do not operate in isolation and so it is unrealistic to as- 
sume that credit events are independent. In reality a whole network of links 
exists between companies in related businesses, industries and markets and the 
impact of individual credit events can ripple through the market as a form of 
contagion. It is thus of fundamental importance when modelling credit, not 
only to understand the drivers of credit risk at an individual company, but also 
the dependence structure between related companies. Whether accounting for 
counterparty risk in the price of a single-name credit derivative or considering 
credit risk in a portfolio context, an understanding of credit dependence is es- 
sential to accurate risk evaluation and pricing. Below it is shown how to do so in 
the case of uncollateralized CDSs on a reference name sold by a risky protection 
seller to a risky protection buyer. 



1.2 Literature overview 



Merton developed the original version of the so-called structural default model, 
which can be viewed as an offshoot of the classical double-entry bookkeeping 
(Merton 1974| ). He postulated that the firm's value a t is driven by a log- 



normal diffusion. The firm, which borrowed a zero-coupon bond with face value 
N and maturity T, defaults at time T if its value ot is less than the bond's 
face value N . Following this pioneering insight, many authors proposed various 
extensions of the basic model 



see, e.g. 



Black and Cox 



Nielsen et al. 1993 , Longstaff and Schwartz 1995 



1976 , Kim ct al. 



Leland and Toft 



1993 



1996 



and Albanese and Chen 2005 among many others. They considered more 
complicated forms of debt and assumed that the default event may be triggered 
continuously up to the debt maturity. One of the main problems with this 
approach is that implied short-term credit spreads are zero given that the default 
time is predictable. In order to avoid this problem and obtain reasonable short- 
time spreads several solutions have been proposed in the literature. It has been 
shown that this can be achieved either by making default barriers curvilinear 
( |Hyer et al.] |1999| , |Avellaneda and Zhu| [200l], |Hull and White| |2001| ), 



jumps into the firm's value dynamics (Zhou 2001b , Hilberink and Rogers 



Lipton 


2002 


and Linetsky 



Sepp 
20rj8fi 



by making default barriers stochastic Finger et al. 2002 , or by incorporating 



or 



2004 , Sepp 2006 , Cariboni and Schoutens 2007 



2002 



Feng- 



Extensions of the structural framework to the two dimensional case have 
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been proposed by Zhou 2001a , Patras 2006 , Valuzis 2008 who considered 



correlated log-normal dynamics for the two firms and derived analytical formulas 
for their joint survival probability using the eigenvalue expansion technique. 
Recently Lipton and Sepp 2009 proposed a novel analytic solution using the 



method of images. In the same paper the authors also propose adding jumps 
to the firm's value processes; this ensures that the default time is no longer 
predictable and solves the problem of zero short-term credit spreads. These 
extensions to two dimensions of the structural model framework have been used 



for the estimation of CVA for CDSs (see for example Lipton and Sepp 2009 



Blanchet-Scaillet and Patras 2011] ). Other approaches, based on reduced form 



modelling, have also been proposed in the literature for this purpose: 


Chen and 


Filipovic 


2003 , 


Leung and Kwok |2005 , 


Brigo and Chourdakis 


2009 , 


Brigo 



and Capponi 2010 and Lipton and Shelton 2011 to mention just a few. 



1.3 Contribution 

The computation of the CVA (DVA) requires studying the joint evolution of the 
assets of the reference name and the protection seller (buyer) in the structural 
framework, provided that the corresponding CDS is viewed from the standpoint 
of the protection buyer (seller) . The simultaneous and consistent calculation of 
the CVA and DVA for a CDS requires the consideration of three-dimensional 
structural models and studying the joint evolution of the assets of the reference 
name, the protection seller and the protection buyer. This task is complex 
both conceptually and technically and, to the best of the authors' knowledge, 
has not been undertaken before. This paper extends the results of |Lipton and| 
Sepp 2009 by considering correlated log-normal dynamics for three firms and 



computing their transitional probability density (the Green's function) for three 
correlated Brownian motions in a positive octant. A semi- analytical expression 
for the Green's function is computed by combining the eigenfunction expansion 
technique with the finite element method. Once the Green's function is known, 
the joint survival probability as well as CVA and DVA corrections for a CDS 
can be computed in a consistent manner. It is worth noting that the proposed 
construction of the Green's function contributes both to mathematical finance 
and to probability theory. 

This paper is organized as follows. Section [2] contains the basic definitions 
necessary for the calculation of credit/debt valuation adjustments. Section [3] 
introduces the structural default model framework. Section [4] shows how to 
price standard single-name credit default swaps in this framework, while section 
[5] extends this calculation to the problem of computing unilateral CVA/DVA for 
standard single-name CDSs. Section [6] contains the main results as it considers 
the three dimensional structural model and obtains a semi-analytical expression 
for the corresponding Green's function. This is then applied to the computation 
of bilateral CVA for a reference-name CDS. The applications of the proposed 
technique to the real market cases are discussed in section[7]where some realistic 
examples of pricing CDSs sold by risky sellers to risky buyers are considered. 
Section |8] gives a brief conclusion. 
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A short version of this paper (Lipton and Savescu 2012| ) has been submitted 
for publication in Risk magazine. 



2 CVA for CDS 

In order to make the paper as self-contained as possible, a brief discussion of a 
standard CDS contract and the corresponding CVA and DVA is presented. By 
entering into such a contract, the protection buyer (PB) agrees to pay a periodic 
coupon c to a protection seller (PS) in exchange for a potential cashflow in the 
event of a default of the reference name (RN) of the swap before the maturity of 
the contract T. The value of a CDS can be naturally decomposed into a coupon 
leg (CL) and a default leg (DL). Let t rn be the default time of the reference 
name, and Rrn its recovery. Then, from the protection buyer's point of view, 
the values of CL and DL are given by: 

CL t = — E [J2 Ti cD (t,Ti) 1 {Ti < t «» } AT| F t ] , (1) 
DLt = E [ (1 - R RN ) D(t, r RN )l {t<T BN <T} | F t ] , (2) 

where are the coupon payment dates and D(t, T) is the price of a zero-coupon 
bond with maturity T. One can simplify the above formulas by denoting by 
CF(t,T) the sum of all discounted contractual cashflows between t and the 
maturity T (both coupon leg and default leg), and writing the value V t of the 
CDS as: V t = E [CF(t, T)\ F t }. 

Assuming now that the protection seller can default but the protection is 
buyer risk free, and denoting by V t the value of the derivative in this case, one 
can represent Vt as follows: 

Vt =E [CF(t, J 1 )l{ r ps >min { T!T flN}} | F t ] 

+ E [l {T P S<min{TiT H« }} [CF(t,r ps ) + D(t,r ps ) (RpsV+ s + V- PS )] | Ft] , 

where t denotes the default time of the protection seller; and, as usual, 
V = ± max (0, ±V). According to the standard market practice, it is assumed 
that if the position is negative in value (to the protection buyer) at the time 
of default of the protection seller, the protection buyer will still be obligated to 
pay in full, while if the position is positive in value they will recover a fraction 
Rps of the value of the position. Due to the fact that V + + V~ = V and 
V t ps = E [CF(t ps ,T)\ F t ps], it can be shown that: 

V t =E [CF(t, T)l{ T PS >mia { TtT RNyy 

+ l{r^ <min{ T,r^} } [CF(t, t ps ) + D(t, t ps )E [CF(t ps , T)\ F t ps] 
- D{t,T PS )(l-Rp S )V+ PS ]\ Ft] . 

Moreover, since D(t,r F ) and CF(t,r p ) are F t ps -measurable, one can write 

CF{t, t ps ) + D(t, r PS )E [CF{t ps , T)\ F t ps] = E [CF(t, T)\ F T rs] . 
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Since t < t ps , it is clear that E [E [ -| J t ps] | F t ] = E[-\F t ] (the tower law). 
Thus 

V t =E[CF(t,T)\F t ] 

+ E [ l {T ps >win{T>rluf}} (CF(t, T) - E [CF(t, T)\F t ps])\ F t ] 
- E [ l {T p S<mln{T , r ^ }}J D(t, t ps ) (1 - R PS ) V+s | F t ] . 

On the other hand one can observe that: 

l{rP S >min{T,r««}} [CF(t, T) — E [CF(t, T) | F t ps]) = 0, 

to obtain: 

V t = V t -E[ l {TPS<min{T , rRN}} D(t, t ps ) (1 - R PS ) V+s | 7i] • 

The term Credit Value Adjustment (CVA) represents the additional cost 
associated with the possibility of the counterparty's default and is defined as 
CVA = y t -t4: 

CVA = (1 - R PS )E [l {T p S < min{ T,^«}}^(^T P5 )F+ s | Ft] . (3) 

Similarly one can consider the case where the protection buyer is risky but 
the protection seller is risk free. The term Debt Valuation Adjustment (DVA) 
represents the additional benefit of one's own default (t pb denotes the default 
time of the protection buyer): 

DVA = (1 - R PB ) E [l {T PB <min{TtT n N}} D(t, t pb )V- pb | F t ] . (4) 

In the current environment, it is no longer reasonable to assume that one 
of the counterparties is risk free. The Basel II documentation makes a clear 
reference to a bilateral counterparty risk, in which both counterparties involved 
in the derivative contract are subject to default risk. This bilateral approach 
introduces much needed symmetry in pricing of a CDS and allows the two 
counterparties to agree on its price (for a detailed discussion on this see for 
example Brigo and Capponi |2010| ). If r denotes the minimum of the two 
default times: r = minlr^, T r±! }, then 

V t =E [CF(t,T)l { T>T} 

+ 1 {t=tPS<t} (CF(t,T PS ) + D(t,T PS )R PS V+ s + D{t lT PS )V; PS ) 
+ 1 {t=t pb <t} (CF(t,r PB ) + D(t,T PB )R PB V T - B + D(t,T PB )V+ B )] . 

In the case where both counterparties are considered risky, bilateral CVA is 
the combination of the two adjustments(CVA and DVA): 

CVA = {l-Rp S )E[l {T Ps <min{T PB^ T n N . T}} D{t 1 T PS )V+ s \ F t ] , (5) 
DVA = (l-Rp B )E[l {TPB<mm{T Ps iT i tNm D(t,T PB )V- pB \ F t ] . (6) 

It should be noted that expressions (§), @ and (|5|, ^ are not identical. 
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3 Structural model framework 



In this section the structural default model for a single name is discussed. For 
simplicity it is assumed that the default and counterparty risk can be hedged, 
so that one can work with the risk neutral pricing measure denoted by Q. It is 
also assumed that cash flows can be discounted with risk-free deterministic rate 
Qt- 

Let at be the firm's asset value. It is assumed that at is driven by the 
following jump-diffusion dynamics under Q (similar to the setup in Lipton and 
Seppl[2059] ): 



da t = (g t - Ct - X t K) a t dt + a t a t dW t 



1) dN u 



(7) 



where Qt is the interest rate, Q is the dividend rate, Wt is a standard Brownian 
motion, <r t is the deterministic volatility, N t is a Poisson process independent of 
Wt, Xt its intensity, j is the jump amplitude, which is a random variable with 
probability density function (PDF) given by u>(j), and k is the jump compen- 
sator: 

r o 

e 3 ui{j)dj - 1. 

3 

Typically, for simplicity, PDFs using one free parameter and negative jumps are 
consider; these jumps may result in random crossings of the default barrier. 

Further, it is assumed that the firm defaults when its value per share becomes 
less than a fraction of its debt per share. In this approach, which is similar to 
that of Finger et al. 2002 and Lipton 2002 , the default barrier of the firm is 



a deterministic function of time given by: 



where 



E t = exp 



h — loEt, 



(8) 



du 



and lo = RLq. Here R is the average recovery of the firm's liabilities (that 
can be estimated from the prices of its bonds and CDS quotes) and L is its 
total debt per share (from the balance sheet as the ratio of the firm's total 
liabilities and the number of common shares outstanding). The convexity term 
|of reflects the fact that the barrier is flat for the logarithm of the asset value 



rather than for the asset value itself (as suggested by Zhou| |2001a and Haworth 



et~al][2008] ) 



Following Stamicar and Finger 2006 , the following approximation of the 



firm's equity price per share s t is used: 



St 



a, 
0. 



t < T 
t > T 



(9) 



where r is the default time. At time t = 0, so is specified by the market price 
of the equity share. Accordingly, the initial asset value is given by do = sq + £o- 
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For simplicity we assume going forward that the volatility is constant in 
time. The solution of the stochastic differential equation ^ can be written as 
a product of a deterministic part and a stochastic exponent: 

at = l Q E t e ax \ 

where the stochastic factor x t is driven by the following dynamics under Q: 

dx t = dW t + 3 -dN u so = - In , (10) 

a a \ < / 

with xo representing the "relative distance" of the asset value from the default 
barrier. In the current formulation, the default event occurs at the first time r 
when x T becomes negative. The default barrier is fixed at zero and the default 
event is determined only by the dynamics of the stochastic driver x%. 



As was emphasized by Zhou 2001b and |Lipton| [2002 , introducing jumps in 



the dynamics of the asset value allows one to calibrate to CDS market spreads 
even for short maturities. In the framework without jumps it is well known 
that the default time is predictable, so that the survival probability is hyper- 
exponcntially flat for very short maturities, and good calibration of distressed 
names in the market impossible. 

The case without jumps however, allows for analytical solutions in some 
cases which are useful for the understanding of the problem, as well as provide 
good benchmark for the more general case with jumps. Besides, CDSs with 
medium and long maturities can be adequately dealt with in the case without 
jumps. Accordingly, this paper is focused on the simplified case without jumps. 

The generalization of the above formulation (especially without jumps) for 
the multi-dimensional case is straightforward. It is assumed that the process 
for the relative distance to default for each of the entities evolves according to 



equation (10), while to corresponding Brownian motions are correlated in the 



usual way, so that d(W^,W^) — pijdt. When jumps are present, stochastic 



drivers can be correlated via a Marshal-Olkin inspired mechanism (see Lipton 
and Sepp|[20fJ9] ). 



4 One-dimensional case 

This section presents the case of the standard single name CDS, where only 
the dynamics of the reference name is modelled, as the protection buyer and 
protection seller are considered non-risky. This case is well known, but discussed 
here for completeness and as a gentle introduction to the subject. 

The process yt measures the relative distance from the default barrier in 
time for the reference name of the CDS. In the simplified case with no jumps it 
has the following dynamics: dy t = dWf, and the starting point y$ > 0. 
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4.1 Pricing problem and Green's function 

The general pricing problem in this framework is given by: 

Vt + \Vyy -QV = 0, (11) 

where the domain is the positive semi-axis: y > Q. Green's function solves the 
forward equation (where r = T — t): 



G T - -Gy> y > — 0, 



with the initial condition G(0,yo,y') — 5 {y' — yo). The solution for this equa- 
tion is well known and given by (using the method of images) : 



G(T,ya,y') = 



1 / (.v'—yp) 2 (y'+vo) 2 



/2tvt 

or by using an integral representation: 

/■OO 



e 2t — e 



G(T,y ,y') = e^^sinfcyo sinky' dk. 

Jo 

Figure [T] shows that the expressions obtained through the two different formu- 
lations coincide. 




Figure 1: Green's function (r = 1 year, y = 1). 



4.2 Survival probability 

We denote by Q(t,T,y) the survival probability to maturity T of the reference 
issuer at time t. This satisfies the equation 

Qt + \Qm = 0, (12) 
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with final condition Q(T,T,y) = 1. Using the Green's function obtained previ- 
ously we can write the analytic formula for the survival probability: 

Q(t,T,y )=[ G(T,y ,y')dy' 



where Af denotes the cumulative normal distribution 



=2Af [ ) - 1, (13) 



4.3 Price of a standard CDS 

In this section we discuss the pricing of a standard CDS on the reference issuer. 
The expression for the coupon leg given in equation ([I]) can be simplified by 
making the assumption that the coupon is paid continuously and using the 
expression in equation ( 13 ) for the survival probability we obtain: 

CL(t,T,y ) = -cA(t,T,y ), 

where A(t, T, y ) denotes the annuity leg and can be written as: 



A{t,T,y ) = / D(t,t')Q(t,t',y )dt'. 



The integral can be computed analytically using integration-by-parts and the 
following indefinite integral (7.4.33 in Abramowitz and Stegun 1964] ): 



i a x * 2 dx = v 



2a 



V2b 



where a ^ 0. The analytical expression for the annuity leg is given by: 



A(t,T,y ) = 



1 



1 - e-^ T ^Q(t, T, y ) - e «°^w(-^= - ^2~g (T - tf) 



yo 



+ ^2g(T-t) 



(14) 



The default leg is given by: 



DL(t, T, yo) = (1 - R RN ) / D(t, t')dQ{t, t' , y ) 



{1-R 



RN) 



1 _ e- e < T -QQ(t,T,y ) - gA(t,T,y ) 



The price of a single-name CDS where both counterparties are considered 
non-risky is: 



V(t, T, y ) =CL(t, T, y ) + DL(t, T, y ) 

= -(c + g(l-R RN ))A(t,T,y ) 



+ (1 - R RN ) 



1 - e 



-g(T-t) 



Q{t,T, y ) 



(15) 
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5 Two dimensional case 



For the two dimensional problem we need to model simultaneously the evolution 
of the asset values for two issuers. Processes Xt and y t measure the relative 
distance from the default barrier in time for each of the two entities considered. 
These processes have the following dynamics: dx t — dWf, dy t — dWf , where 
the Brownian motions W x and W y are correlated with correlation p xy , \p xy \ < 1. 



5.1 Pricing problem 

The general pricing equation in this framework is given by: 

V t + ~V XX + ^V yy + p xy V xy -qV = 0. (16) 

We consider the following function U(t, T, x, y) — e e ^ T ~^V{t, T, x, y) and apply 
a change of variables that allows us to eliminate the cross derivative and killing 
term: 

a(x,y) = x 

y) = - ^ (p xy x - y) , 

Pxy 

where we have used the notation: p xy = yjl — p xy . This leads to the following 
simplified version of the pricing equation: 

U t + \u aa + hjpp = 0. (18) 

Along with the change of variables, the domain this has to be solved in has 
changed from the positive quadrant to the interior of an angle (see figure [2J). 
This angle is characterized by cos(iu) = —p xy , so if p xy > 0, the angle is obtuse. 

In order to take advantage of the symmetry of the domain, we make a second 
change of variables and convert to polar coordinates: 



I a = — r sva\(p — w) I 

\P=r cos(</? -w) y = w + arctan 

The final form of the pricing equation becomes: 



(19) 



U t + \ \ U rr + -U r + ^U VV ]=U. ( 21) ) 

5.2 Green's function 

We concentrate now on calculating Green's function by solving the forward 
equation: 

G T - - \ G r i r > + —,G T : + ^G^'v'j = 0, (21) 
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Figure 2: The new domain in which the PDE has to be solved in after the 
change of variables. 



with initial condition: 

G(0,rV) = l<Kr' - r )%' - <po), 

and zero boundary conditions: 

G(r,r',0) = 0, G[r, r', w) = 0, G(t, 0, tp') = 0, G(r, r', if') > 0. 



The polar coordinates (ro, tpo) of the source are given by: 



''0 



y/xl ~ 2p xy x y + yl 
— i 

Pxy 

ipo = arccos (—p xy ) + arctan 



Pxy X 
2/0 - PxyXQ 



Two possible methods can be applied in order to obtain the solution for 
Green's function: the eigenvalue expansion method and the method of images. 
The solution using the first method is well known and has been first introduced 
in He et al. 1998 , Lipton 2001 , Zhou 2001a . We give a brief outline in section 



5.2.1 A solution through the method of images was announced by Lipton in 
2008 at a SIAM meeting, and briefly discussed in Lipton and Sepp 2009). We 



give in section |5.2.2| a detailed presentation on how to obtain Green's function 
through this method. 



5.2.1 Eigenvalues expansion method 

In this section we aim at giving a solution for Green's function through the 
eigenvalues expansion method. This is a well known method for solving Green's 



equation and has been extensively studied in the literature (He et al. 1998 
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Lipton| |2001| , |Zhou| |2001a] , |Patras| |2006| , |Valuzis| |2008| ). We give a brief 
outline of the methodology here as it is instructive and a starting point for the 
new methodology we develop in section |6.2| for the three dimensional case. 
The first step is to apply the separation of variables technique: 

G{T,r',ip') = g(T,r')f(ip'), 

where the zero boundary conditions for the Green's function now apply to func- 
tion /: /(0) = and f(zu) = 0, while for the function g we have the initial 
condition: g(0, r') — £-6 (r' — Tq) and boundary conditions g(r, 0) = and 
g(T,r')- >(). 



By substituting back in equation (21 ) we can rewrite the equation such that 



the left hand side depends only on r and r' and the right hand side depends 
only on ip'. Hence both sides are equal to some constant value C and we have: 

C 



1 



9t 



1 



., i 9r'r' T . 9r' 

2 V r 



--Cf. 



It is well known that we necessarily have C < and hence we make the 
notation C = —A 2 . Imposing the boundary conditions for function / we obtain 
that A = ^ for positive integers n, and the solution is given by / (ip 1 ) = 



. We now proceed to solving the PDE for <?(r, r'): 



9t 



9r 



1 

-,9V 



A 2 



(22) 



with the corresponding initial and boundary conditions. We claim that the 
solution is given by: 



-Ia 



r r 

T 



where Ia(0 is the modified Bessel function of the first kind and satisfies the 
following equation: 

,d 2 I dl 



(C 2 +A 2 ) =0. 



(23) 



One can verify that this is indeed the case by computing the relevant derivatives 



of g and substituting back in equation (22 ). We can also verify that the function 



satisfies the initial condition. For this we use the asymptotic approximation for 
the modified Bessel function in the limit where £ 3> A: 

pi 

Ia(0 



/27T£' 



and we obtain: 



g(i-,r') 



> —5 (r' - r ) . 
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Similarly we can show that function g also satisfies the boundary conditions at 
r' — » and r' — > oo. 

Having solved separately the equations obtained when applying the method 
of separation of variables we can now write the solution for the Green's function: 



-/2 , 2 
r +r 



G(t, r ,r', ip ,ip') 



r r \ I nirtp 
sin 

T J \ W 



To simplify the equations we use the following notation: v n = — . The 
coefficients C n can be computed by imposing the initial condition for Green's 
function and we obtain that 

oo 

^ C n sin (i> n tp') =5(ip' - ip ) . 
n=l 

We multiply by sin (y m ip') and integrate from to w 1 and we have for the 
coefficients the following expression: C n = ^sm(v n tf'). The final formula for 
Green's function in the domain shown in figure [2] is: 



G(t, r ,r',ip ,(p') 



2e~ 



n=l 



r r a 



sin (v n (p') sin {v n y Q ). (24) 



Figure [3] shows the two dimensional Green's function for sample values for 
the input parameters. 



Green's Function 




Figure 3: Green's function (x a = 0.1, y a = 0.1, a x = 10%, a y 
70%, T = 1 year). 



10%, p. 



13 



5.2.2 Method of images 

In this section we aim to give a solution for Green's function through the method 
of images. This has been announced in |Lipton and Sepp 2009 and we give here 
a detailed presentation on how to apply this method in our case. 

We first need to find the solution to equation (21) with the same initial 
condition but with non-periodic boundary conditions: 

H(t, 0,<f/) = 0, H( r y,<p') y0, H(T,r', V ) ► 0. 

r'— >oo W\— 

We perform the Fourier transform in if and denote by H (r, r' , v) the shifted 
Fourier transform of H (t, r',tp'): 

/CO 
H{t, r',Lp')e- l ^'dv'. 
-co 

We obtain the following problem for H (r, r' , v): 

H T - - (h^ + —H r , - ^H^j = 0, 
with boundary conditions i?(r, 0, </?') = 0, H(r,r' ,ip') > and the initial 

r' — fco 

condition: H (0,r' \v) — i#(r' — r ). We recognize that this problem is the 
same as in equation ( 22 ) . We have shown there that the solution to this equation 
is given by the following expression 

- r +r O / I 

til i \ e 2T t I r r ° 
H{t, r ,v) = I\ v \ I — 

In order to obtain H (r, r' , ip') for the problem with non-periodic boundary 
conditions we use the inverse Fourier transform and represent it as an integral 

H(T,r S,<Po,</) = —^Tr J *M r-f)^-^dv 

— oo 

3 -(r' 2 +r 2 )/2. ; , r(i 

I v ( — — J cos [y (ip' - ipo)) dv. 



ii 



We observe that the integrals depend only on the difference if' — ipo, which 
we denote by if). To simplify this formula we use the following integral repre- 
sentation of the modified Bessel function l v (z) for nonnegative i/, z: 

7T OO 

I v {z) = i / e zcos0 cos{vd)de- [ e- zcosh ^dC 

7T J 7T J 
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Accordingly, we can write H as follows 

e -0' 2 +^)/2T 

H (t, r , r', V) = 



.4 



r r 

T 



B 



t r 

T 



where the functions A(z) and B(z) are defined as: 



A{z) 



— / cos (v6) cos (v4>) dv 

TT 



z cos 6* 



oo oo 



B (z) = J — J e sin (z/7r) cos ^ 



o L o 



d0, 



3 — 2 cosh £ 



The inner integrals with respect to z/ can be easily calculated and we obtain 

1 



4( Z )= 2 e zras *l { ^< f) , 



[C 2 +7T 2 - e-* cosh(; 



d(. 



C 2 + (tt + ^)"J C 2 + (tt - V) 
Finally, we obtain the following expression for the Green's function 

e -(r' 2 -2rV cos(V)+rg)/2r 

H(T,r ,r',ip) = — J [-7r,7r] ("0) 



1 

7rr 



f^2 _|_ ^2 _ ^2j e -(r' 2 +2rV coshC+r 2 )/2r 



C 2 + (tt + ?A) 2 C 2 + (tt - V) 2 



d< 



=i?i (r, r 0) r', VO - H 2 (r, r , r', V) . 



Integral B{z) is discontinuous at ip = ±n (as it can be seen in figure [4]). 
However, even if H changes its form along the lines tp — ±7r as a consequence of 
these discontinuities, it can easily be verified that it is smooth and well-behaved 
(see figure [5]). 

We can transform H lt2 (r, r' , ro, ip) as follows 



Hi (r, r ,r', *) = e ' ( ?"° )/2r(5+ : S - ) e (^/-) «W , 



#2 (t, r ,r',ip) 



e 2t 



2ttt 



27T 2 7 



C 2 +0+«/<) 1 c 2 +(^-iW 



-(r'ro/r) coshC^/- 



e 2t 



27r 2 r 



S+e" 1 ^ cosh((7T+VX) + s _ e -^a cosh(( 7 r-^.)C) 

C 2 + i 
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phi-phi' 



Figure 4: Discontinuities of the function B(z) at ip = ±7r for z = 1. 
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where s± = sign (jr ± ip) . We want to rewrite the above expressions in a more 
compact form. To this end we introduce the following function / (p,q), where 
p > 0, — oo < q < oo: 

1 f°° p-P( cosh ( 2 9C)-cos((j)) 

and its extension h(p,q) is defined as follows: 

h (p, 9) = 2 ( p > ^ + ?) + s -/ (P> i" - «)] • 

Then we can represent H (t, ro, r', "0) in the following form, which can be viewed 
as a direct generalization of the one dimensional case: 

H(T,r y,ip) = + ° 2 - W> ° h \ 

Z7TT \ T 

Now that we have obtained the solution for the problem with non-periodic 
boundary conditions, we can go back to our problem of interest which requires 
us to solve equation (21 1 in an angle, where <<£>'< tu. For this problem we 
can represent the fundamental solution in the form 

00 

(t, r , r', ip , ip') = (r, r , r', ip + 2nru, <p') 

71= — OO 
OO 

- ^ f H(T,r ,r f ,-<po + 2nw,<f/). (25) 

n— — oo 

Indeed, it is clear that these sums converge, every term solves the parabolic 
equation and only one term has a pole inside the angle. After obvious rear- 
rangements, we can write: 

00 

H m {r, r , r\ ip , 0) = ^ [H (r, r , r', ip + 2nw, 0) - H (r, r , r', -ip - 2nru, 0)] = 0, 



H m (T,r o y, ip Q ,w) = ^2 i H (t^o,^,^ + 2nw, w) 

n— — oo 

—H (r, ro, r , — ipo + 2vo — 2nw, m)] = 0, 

by symmetry. As expected, results using the representation given in equation 
(25) coincide with those obtained using representation (24 1, obtained through 
the eigenvalue expansion method. 
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5.3 Joint survival probability 

We denote by Q(t, T, x, y) the joint survival probability of issuers x and y to a 
fixed maturity T. This solves the following pricing equation 



with final condition Q(T,T, x,y) — 1 and boundary conditions Q(t,T, x,0) = 
and Q(t, T, 0, y) = 0. After applying the change of variables described in section 
|5.1[ we obtain the following PDE: 



+ 



1 



r 



Q 



2 ^(pip 



= 0, 



with final condition Q(T, T, r, ip) = 1, and boundary conditions Q(t, T, r, 0) = 
and Q(t,T,r,w) = 0. We use the expression for Green's function obtained 
through the eigenvalue expansion method and we obtain for the survival prob- 
ability (r = T-t): 



Q(r,r ,ip ) 



2r'e~ 



o jo 

oo 



VJT 



— E ( ! T £i ) s Mvn<p') sw,(i/ n (p )dip' dr' 



E 

k=0 

oo 

E 

fc=0 



(2k + 1) 7TT 



sin(z/ 2 fc + i</>o) I re ^ I„ 2k+1 ' I 



V T / 



4 sin(^ 2fc+1 y ) / rg V2fc+ir (l+ 

(2fc + i)7r ^2ry 2 r(i+^ 2fc+1 ) 



(26) 



where i-Fi denotes the confluent hypergeometric function. This last expres- 
sion allows for a generalization to the three dimensional case, which we discuss 
later. For the two dimensional case this can be simplified further (for details 
Iyengarj [1985] or |Metzler| [2010] ): 



see 



s 2r e 4t v-^ sin (^k+ifo) 
Q{T,r , tpo) = k — 2^ 



k=0 



2k 



(27) 

Figure [6] shows the joint survival probability for two issuers for a range of starting 
point values and two sample correlations. 



5.4 Application to the CVA computation 

We associate the process x t with the protection seller and the process y t with 
the reference name issuer of a CDS. The protection buyer will be considered 
non-risky in this case. The pricing equation for computing the CVA is given by: 
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Joint survival probabilty: corr-75% 



■ 80.0%-100.0% 

□ 60.0%-80.0% 

□ 40.0%-60.0% 

■ 20.0%-40.0% 

□ O.u%-20.0% 




^-^■COCOOJt-^-OO 

odoodoooo 



Inital value x 

(a) p xy = 75% 



Joint survival probabilty: corr--75% 




Inital value x 



(b) p xy = -75% 



Figure 6: Joint survival probability for sample correlations and starting points 
values (a x — 15%, a v — 15%, T — 1 year). 



V t + ^V XX + -V yy + Pxy V xy -qV = 0, (28) 

with the final condition V(T, T, x, y) = and the boundary conditions depend- 
ing on the payoff. In the case of the CVA calculation these are: 

• If the credit referenced by the CDS contract defaults first: since the pro- 
tection seller has not defaulted it will be able to honour the payment and 
hence we have V(t, T, x, 0) = 0. 

• If the protection seller defaults first: it will no longer be able to honour 
its payments and hence the shortfall for the protection buyer will be a 
fraction of the outstanding present value of the single name swap: 

V(t, T, 0, y) = (1 - R PS ) V CDS (t, T, y)+ . 

• If the protection seller is risk free there is no shortfall: V(t, T, oo, y) = 0. 

• If the CDS reference name is virtually risk-free we do not care what hap- 
pens to the protection seller: V(t, T, x, oo) = 0. 

After applying the function and first variable change as in section |5.1[ the 
pricing equation becomes: 

U t + ^U aa + hj w = 0, 
with final condition U(T, T, a, ft) = and boundary conditions: 

U[t,T,a,-^-a ) =0, U(t, T, oo, j3) = 0, U(t, T, a, oo) = 0, 

V Pxy J 

U{t,T,Q,0) = e e(T - t} (1 - Rps)V cds (t,T,p^/3) + . 
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Applying next the second change of variables given in equation ( 19 ), we have 
the following pricing equation: 

U t + \ \Urr + \u r + ^U v A = 0, (29) 

with the final condition: U(T, T, r,<p) — and boundary conditions: 

U(t,T,0,tp) = 0, U(t,T,oo,tp) = 0, [7 (i, T, r, 0) = 0, 
U(t, T, r, w) = e^ T -*) (1 - i? PS ) ^ CDS (t, T, ^r) + . 

In order to obtain the solution [/ that satisfies the pricing equation (291, we 



start from the following identity: 

T oo w _ 

C/f + ^ (u rr + ^U r + ^U vv ^j G(t' - t, r, ip)rd<pdrdt' = 0, 

too" 

and perform a series of integration by parts. We then use the boundary condi- 
tions, the initial condition for Green's function, and final condition for U, along 



with the fact that Green's function solves the forward equation (21), and we 
obtain the final solution for our problem: 

U(t,T,r Q ,ip ) = 

Too 

' 1 [G tp (t'-t,r,0)U(t',T,r,0)-G v (t'-t,r,w)U(t',T,r,zu)]-drdt'. 



2 



t o 



We specialize this expression for the boundary conditions we have for the 
CVA problem and we get: 

Too 

V OVA (t,T,r , l p ) = -' i —^J jD(t,t')G v (t'-t iri uj)V CDS (t\T,p xy r) Urdt' . 

t o 

(30) 



6 Three dimensional case 

For the three dimensional problem we need to model the dynamics of the asset 
values of the reference name, protection seller and protection buyer simultane- 
ously. Processes Xt, yt and z t measure the relative distance from the default 
barrier in time for each of the three entities considered. These processes have the 
following dynamics: dxt — dWf, dy t — dW^ , dz t — dWf, where we correlate 
the Brownian motions with correlations p xyi p xz , p yz . 
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6.1 Pricing problem 

The general pricing problem in the octant: 

V t + -V xx + ^Vyy + ~V ZZ + P x yV x y + PxzV xz + PyzVyz ~ QV = 0. (31) 

We consider the following function U(t,x,y,z) — e e ^ T_t ^(i, x, y), and intro- 
duce a change of variables that allows us to eliminate the cross derivatives: 

' a(x, y, z) =x 
P(x,y,z) {-Pxyx + y) 

Pxy K 6 ^) 

7(ar, y, z) =^— [{p xv Pyz - Pxz) x + (p xy p xz - p yz ) y + p 2 xy z] , 

PxyX 



where we use the notation x — y 1 — Pxy ~ Pxz — Pyz + ^PxypxzPyz- In order 
for the change of variables to be valid we consider \p xy \ < 1 and p xy , p xz and 
2 p 2 

cy rxz "yz 



p yz such that 1 — p\ y — p xz — p yz + 2p X yp xz Pyz > 0. The equation we need to 



solve simplifies to: 

u t + -u aa + X -Ufi (i + i[/ 77 = 0. 

With the change of variables, we have also changed the domain in which we 
need to solve the pricing problem. The original domain was the volume bounded 
by the planes x = 0, y = and z = 0. This now changes to the volume bounded 
by the planes: a = 0, (a,-^a, 7 ) and (a,/3, ?f (-p xz a + P ' vP £~ Pv ~ p) ) i 
we denote by IT, n 2 and n 3 respectively the three planes. We denote by e*3 the 
versor corresponding to the edge IT n II2 , by e% the versor corresponding to the 
edge ITi n II3, and by e*i the versor corresponding to the edge 112 0113: 

§8 =(0,0,1), 

-> / n X Pyz ~ PxzP X y 

e-2 = U, - — — 



ei 



PxyPxz PxyPxz 
X PxyX Pxz ~ PyzPxy 



\Pyz Pxy Pyz Pxy Pyz 

The domain of interest has changed to the hull spanned by these vectors: 

v = u)ie[ + w 2 e"2 +w 3 e3, Wj > 0. 

In order to take advantage of the symmetry of the problem we perform a 
second change of variables to spherical coordinate^] the axis a — and = 
is given by 9 = 0; the axis a = and 7 = is given by tp = and = ir/2. 



a =r sin # sin 95 
/? =r sin 9 cos 
7 =r cos 



r =^a 2 + (3 2 + 7 2 
= arccos ^— j 

(/3 =arctan 
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Y 




Figure 7: Domain after the change in coordinates for p xy — 20%, p xz = 0%, 
p yz = 30% 

In order to obtain the range of possible values for ip for the domain of interest, 
we project e[ and onto the (a, /3) plane and obtain the following normalized 
vectors: 

H 2 = (0, 1) , 

H l = {Pxy -Pxy) ■ 

The range of values for ip is therefore given by: < (p < w, where 
w = arccos (— p xy )- As can be observed in figure [TJ the possible range of 
values for depends on ip: < 9 < 9(y>). 

In order to calculate (ip) we first consider a vector on the boundary of the 
domain (in the II3 plane): 

X = - (uip yz ei+p xz e 2 ) , 
X 

where u> > (the constants are added for convenience in the calculations). Using 
the formulas for e{ and €2, we have: 

^ / 1 - PxyU U) (p xz - PyzPxy) + Pyz ~ PxzPxy 
A — CJ, — 



U)p xy 1 - p xy UJ 



Pxy XPxy 

The projection of this vector onto the (a, j3) plane is the following (normalized) 
vector: 



a/1 - 2p xy u: + w 2 ' y/l - 2p xy oj + i 



1 Notice that the change to spherical coordinates is not the classical one since tp = and 
8 = 7r/2 denote the /? axis rather than the a one. This is done for convenience such that the 
range of possible values for ip is between and a maximum value. 
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and we obtain the angle ip as a function of uj: 

1 - p xy u> 



ip (w) = arccos 



2p xy (jJ + U! 2 



(33) 



It is easy to verify that this parametric form for ip has the right bounds: tp(Q) = 
and tp (w) > vj. In order to obtain a parametric form for 9 we compute the 

UJ — >oo 

length of the vector X which is given by: 



X = 



yl - PL - 2w (Pxy - PxzPyz) + OJ 2 (1 - P 2 yz ) 



and we obtain: 

8 (w) = arccos 

In particular we have 



Pyz PxzPxy ~t~ ^ (,Pxz Pyzpxy) 
\J~Pxy (Plz ~ 2w (Pxy ~ A^Pyz) + W 2 p^ 2 ) 



(34) 



9 (0) = arccos 
6H - 



> arccos 



Pyz PxzPxy 
PxyPxz 

Pxz PyzPxy 
PxyPyz 



Formulas (33 ) and (34) give a parametric characterization of the boundary of 



the domain which will prove very useful going forward. In the domain described 
above, the final form of the pricing equation is given in equation (35): 



u t 



i 



i 92 

r dr' 



(rU) 



1 



1 



r 2 Vsin 2 



:U, 



1 d 
s^9d9 



0, 



(35) 



with appropriate boundary conditions depending on the payoff we are interested 
in. 



6.2 Green's function 

We now concentrate on solving the forward equation for Green's function in 
spherical coordinates: 



1 



r or 2 - r A \ sm 



1 Id 

9 G ^ + ^de> {s ™ e ' G9,) 



0, (36) 



G(0,r',if 



' .J a'\ 



r 2 sin 6 



6(r' -r )5(p' -^ )5(9' -9 ). 



G (r, r', 0, 0') = G(r, r' , w, 9') = G (r, r', ip' , 0) = 0, 

G(r, r', tp', 6(^0) - G(r, 0, ip' , 9') = 0, G(r, r', y/, 0') ► 0. 



23 



We aim to build a solution for Green's function through the eigenvalues ex- 
pansion method. The first step is to apply the separation of variables technique: 

G{T,r'^',e') = g{ry)n<p',n (37) 



By substituting (37 1 into (36), we obtain an equation where the left hand side 
depends only on r and r' and the right hand side depends only on ip' and 9' , 
and hence both sides are equal to some constant value C, which is necessarily 
negative. We use the notation C = —A 2 , and obtain the following equations for 
functions g and VP: 

1 f 1 d 2 , . . A 2 N 
( r 9) ~ -j^9 



1 Id 

^-aZ7*«>V + "^7 7m7 (sin0'*eO = -A 2 *, 
sin 0' sm 0' ay' 

For function g(r, r') we have the initial condition g(0, r') = \8 (r' — r ) and 

boundary conditions g(r, 0) = and g(T,r') > 0, while for function 9 we 

have zero boundary conditions: 

9(0,0) = 0, v(za,6') = o, *(^,o) = o, *(^,e(^)) = o. 

6.2.1 Radial part 

To solve the PDE for g(r, r'), we introduce a new function h(r, r') — y/r'gfr, r'). 
The equation that h satisfies is: 

If 1 A 2 + 1 \ 
h T = -[h r , r , + -h r ( 38 ) 

with the initial condition h(0, r') = rg ^— 8(r' — ro). We observe that this is 
a similar equation to equation ( |22[ ), which was solved for the two dimensional 
case. Similarly to that, the solution for this equation is given by 



h(r, r 1 ) = - — — / 



rrQ 

which yields the following solution for g: 



6.2.2 Angular part 



In order to obtain Green's function for the desired problem, we also need to 
solve the two domensional PDE for 9(>p', 9'): 

V^ + -474( sin ^') = -A 2 *, (39) 



sin 2 ^ 1 ^ ' sin 0' 90 
*(0,6»') = 0, V(m,6') = 0, Q) = 0, 9(ip', Q(<p')) = 
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The eigenvalue problem given in equation ( 39 ) but considered on the surface 



of the whole sphere is a well known problem. It has been shown, in Courant and 



Hilbert 2008 for example, that this problem has a countably infinite sequence 
of positive eigenvalues < A x < A 2 < as well as a corresponding sequence 
of linearly independent eigenf unctions. The solutions are obtained using the 
separation of variables technique and are known as the spherical harmonics. 

However, in our case, a further separation of variables is not possible because 
of the particular shape of the domain. The two dimensional spherical surface 
inside the red line shown in figure [7] can be mapped directly onto the (ip',6') 
plane. This is done in a similar way to the method used by cartographers to 
map the Earth's surface using Mercator's projection. The southern boundary of 



the domain is mapped into a continuous curve parametrised by equations ( 33 1 
and ( [34| . The boundary at 6' = is degenerate as it corresponds to the north 
pole on the sphere. 

Figure [8] shows the domain (denoted hereafter by it) projected onto the 

show 



11 



(ip',6') plan when all correlation values are set to 0. Figures [9 10 and 
the oriented domain for sample positive correlation values, while figures 12 and 
13 show the domain for sample negative correlation values. 
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Figure 8: p xy 
0, p yz = 



0, p x 



Figure 9: p xy 
0.5, p yz = 0.5 



0.5, p xz 
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Figure 10: p xy 
0.5, p yz = 0.3 



= 0.8, 



Figure 11: p xy 
0.05, Pyz = 0.6 



= 0.8, 




Figure 12: p xy = 0.2, p xz = Figure 13: p xy = 0.8, p xz = 

-0.1, Pyz = -0.6 -0.65, Pyz = -0.45 

Given the varied forms that the boundary of the domain can take, as well 
as the fact that it is a curved boundary, we construct the solution to this 2D 
PDE using a finite element method. 

We note that '5 should satisfy the same regularity conditions over the whole 



surface of the domain. As noted in Courant and Hilbert 2008 , the operator in 



equation ( 39 1 is invariant under rotations of the coordinate system. Therefore 
the singularity at the point 8' = is a consequence of the particular choice of 
the coordinate systemj^] 

2 Since for our domain we have 8' < it (the south pole can only be reached when p^. y = 1 
which has already been excluded in order for the change of variables | |32| to be valid), we 
could make an infinitesimal rotation of the syst em o f coordinates such that through our whole 
domain < e < 9' < ir holds. Since equation ( |39| is invariant under rotations, we have the 
same eigenvalue problem but on a domain that does not contain the singularity and hence 
our solution will have all the required regularity properties. Since e is arbitrarily small, the 
change in the domains [8| to [13| is not noticeable. 
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In order to obtain the variational formulation (or weak formulation) of the 



spectral problem given in (39 1, we use a test function and integrate over 
the whole domain. The test function ^f' belongs to the same space as in 
particular it is also on the boundary of the domain. Using integration by 
parts and Green's theorem, along with the fact that the test function is null 
over the border of the domain, we obtain the weak formulation: 

/ -*„,/*', dfi+ / sm6'^ e ,%,dn = A 2 [ Wsin.e'dn. (40) 

J n sm0' v Jn Jn 

To obtain the solution through the finite element method, we start by con- 



structing a triangular mesh for the domain (section 6.2.3 describes in some 
detail how the mesh is built). The space in which we are searching for the so- 
lutions is replaced by a finite dimensional space. The dimension of this space is 
given by the number of free points in the mesh, denoted by n (the number of 
vertices of all triangles in the mesh excluding those that are on the boundary 
of the domain). The finer the mesh is, the higher the dimension of this space, 
and the better the approximation of the solution is. We denote by ($i) 1 < i <„ 
the basis for this space. 

We consider linear basis functions on each triangle, and given any triangle 
T of the mesh, there are only three basis functions that are non-zero on T. Wc 
denote by (cp' lt (<p' 2 , 9' 2 ) and (ip' 3 , 6' 3 ) the vertices of triangle T and by $i, $2 
and <j> 3 the corresponding non-zero basis functions. These are defined by: 



*< (^i) 



1, i= j 
0, %+3 



and can be represented by $>i(<p', 9') = a, + bi<p' + Cj0, (cp', 9') e T, i = 1,2, 3. 
The coefficients a, , and Ci can be found by solving a system of 3 x 3 equations. 

The solution for our problem can be associated with a vector in R n and can 
be written as: ^ (<//, 6') = Y^i=i (v 5 '; The weak formulation given in 



( 40 ) is then approximated by the linear system: 

A"<5 = A 2 M*, 

where we denote by K = (Kij) 1<i - <n the stiffness matrix, and by M 
(-&%)i<jj< n the mass matrix: 



Kij = [ (AV<f>j) ■ V^dft, 
Jn 

I sine' Q^jdn, 

Jn 



with the matrix A given by A = si " e ' ././). Each of the integrals in- 
J V sin0 J 

volved in the computation of the elements of matrices K and M can be rewritten 
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as a sum of integrals over the triangles where the basis functions are non-zero: 
Kij =Y^J CAVSj) • V$idT k 

Since the basis functions are linear over the triangles where they are non-zero, 
the derivatives are constant, and hence the computation of the elements of K 
comes down to the computation of integrals J T ^jydT k , J T smO'dT^ over the 
triangles in the mesh. This can be done by the standard "one-point" quadrature 
rule, for example: 

f f{<ff,ff)dT k =Ana Th f 

where ((/?, #) is the centroid of triangle Tj. (higher precision quadrature rules can 
be used as well). The elements of matrix M can be computed in a similar way, 
and we can now solve the linear system associated with the weak formulation 
of our problem. 

To solve this linear system we first do a Cholesky decomposition of matrix 
M (note that the matrix M is symmetric): M = M.M. T and we have: 

M^K^> = A 2 M T y. 

We introduce matrix C defined by C = Mr x K (A4 _1 ) , which is also symmet- 
ric, and the system can be rewritten as: 

C(M T V) =A 2 (M T ^). 

We compute the eigenvalues and eigenvectors for this problem, and the eigen- 
vectors for the original problem can be computed as \-M T ) E where E is an 
eigenvector of the modified problem. Sample results for this eigenvalue problem 
in our particular domain are discussed in section |6.2.4| 



6.2.3 Constructing the grid 

We give here a brief description of the methodology used to construct triangular 
meshes on the domain of interest. The algorithm used is iterative. The nodes 
of the mesh are adjusted at each iteration based on the current element sizes 
according to the ideas presented in |Persson| [2005] . The Delaunay triangulation 
algorithm is then used to adjust the topology (decide the edges) at each iteration. 
For the Delaunay triangulation we use a divide and conquer algorithm along with 



the quad-edge data structure described in detail in|Guibas and Stolfi 1985 



Figures [14] and [15] show the uniform meshes obtained with this method for 
two sample sets of correlations. 
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(a) First iteration mesh 



(b) Mesh after 100 iterations 



Figure 14: Uniform mesh for the domain obtained for p xy — 0%, p xz = 0%, 
p yz = 0%. The mesh is constructed using 1500 points. 




0.5 1 1.5 2 2.5 0.5 1 1.5 2 2.5 



(a) First iteration mesh (b) Mesh after 100 iterations 

Figure 15: Uniform mesh for the domain obtained for p xy = 80%, p xz = 20%, 
Pyz — 50%. The mesh is constructed using 1800 points. 

However, there are cases where it is advantageous to have different sized ele- 
ments in different regions: where the geometry is more complex or the problem 
requires more accuracy (for example close to a singularity such that the global 
accuracy of the solution is good). In order to create adaptive meshes for our 
domain, the desired edge length distribution over the domain can be specified 
(this does not have to equal the actual size, but it rather gives the relative 
distribution over the domain). 

Algorithm [T] gives a brief description of the method used to build adaptive 



2') 



Algorithm 1 Algorithm for constructive an adaptive mesh 



Require: X%, X2, Y\, Y-x - bounding box of the domain 

Require: d(x, y) - distance function to the closest boundary (negative inside 
the domain) 

Require: h(x,y) - element size function (gives the relative element size distri- 
bution over the domain) 
1: Build a mesh with equally spaced points for the bounding box of the domain 



Remove points outside the domain 

Rejection method: reject points inside the domain with probabilities pro- 
portional to l/h(x,y) 2 
while i < MAXITER do 

Delaunay triangulation using the divide and conquer algorithm described 



in detail in Guibas and Stolfi 1985 



Assemble triangles obtained through the Delaunay procedure 
for each triangle do 
Compute centroid 

If centroid outside the domain: remove triangle from list 
end for 

Move mesh points based on current edge lengths using ideas described in 



Persson 2005 



12: Bring points that have moved outside of the domain back to the boundary 



13: i = i + l 
14: end while 
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triangular meshes. The fixed number of iterations can be replaced by a condition 
on the largest move of a point in the mesh in the previous iteration. 

To obtain the uniform meshes shown in figures [14] and |15[ the element size 
function is constant over the domain. This means that step 3 in algorithm [T] 
does not reject any points. Figures [T6| and [T7| show examples of meshes obtained 
for different values of the correlations and a non-uniform element size function. 
The meshes are finer close to some or all four of the boundaries. 

In each case a mesh similar to the ones used as starting point for the uniform 
case is constructed first (by performing steps 1 and 2 in algorithm [T]) . Then 
the rejection method eliminates points in the regions where we do not need as 
much precision. The Delaunay triangulation of the remaining points is used 
as the starting mesh for the iterative process (steps 4-14), and is denoted in 
the graphs as the "first iteration" mesh. The figures also show the final mesh, 
obtained after 100 iterations. Figure [18] shows a similar example, when two of 
the pairwise correlations are negative. 




Figure 16: Adaptive mesh for the domain obtained for p xy — 80%, p xz = 50%, 
pyz = 50%. The mesh is constructed using 1500 points. The mesh is finer near 
two of the boundaries. 
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(a) First iteration mesh 



(b) Mesh after 100 iterations 



Figure 17: Adaptive mesh for the domain obtained for p xy = 80%, p xz = 50%, 
Pyz — 30%. The mesh is constructed using 1500 points and is finer as we get 
closer to the boundaries. 
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(a) First iteration mesh (b) Mesh after 100 iterations 



Figure 18: Adaptive mesh for the domain obtained for p xy = 20%, p xz = —10%, 
p yz = —60%. The mesh is constructed using 1600 points and is finer as we get 
closer to the boundaries. 

6.2.4 Eigenvectors 

Once the mesh is constructed, we solve the eigenvalue problem in matrix form 
and we obtain the eigenvalues and corresponding eigenvectors. Figure [T9| shows 
the case where all correlations are 0. Figure [20] shows sample eigenvectors 



for a domain where all three correlations are positive, while figure 21 shows a 
case where two of the correlations are negative. Even though the shape of the 
domain varies significantly between the different examples, we observe the same 
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patterns, with an increasing number of modes for higher order eigenvectors. 
Note also that for the first eigenvectors the modes are better defined than for 
the higher order ones. 
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0.2 0.4 0.6 0.8 1 1.2 1.4 

(a) Eigenvector 1: Af = 12.0 




0.2 0.4 0.6 0.8 1 1.2 1.4 

(c) Eigenvector 3: A§ = 30.2 




0.2 0.4 0.6 0.8 1 1.2 1A 

(e) Eigenvector 10: Af = 92.4 
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0.2 0.4 0.6 0.8 1 1.2 1.4 

(b) Eigenvector 2: A| = 30.2 
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-15 
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(d) Eigenvector 6: A| = 56.8 




0.2 0.4 0.6 0.8 1 1.2 1.4 

(f) Eigenvector 20: A| = 189.2 



Figure 19: Eigenvectors and corresponding eigenvalues for the domain obtained 
for p xy = 0%, p xz = 0%, p yz = 0%. The mesh is constructed using 1500 points. 
The mesh is finer as we get closer to the boundaries. 
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0.5 1 1.5 2 2.5 

(a) Eigenvector 1: = 5.2 




(c) Eigenvector 3: A| = 16.3 




0.5 1 1.5 2 2.5 

(e) Eigenvector 8: A| = 39.1 



0.5 1 1.5 2 2.5 

(b) Eigenvector 2: A? = 11.8 




0.5 1 1.5 2 2.5 

(d) Eigenvector 4: A| = 21.3 




0.5 1 1.5 2 2.5 

(f) Eigenvector 30: A§ = 140.0 



Figure 20: Eigenvectors and corresponding eigenvalues for the domain obtained 
for p xy = 80%, p xz = 20%, p yz — 50%. The mesh is constructed using 1800 



points and is shown in figure 15 
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0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 

(a) Eigenvector 1: Af = 21.5 



0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 

(b) Eigenvector 2: A| = 42.2 
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(c) Eigenvector 3: A§ = 63.8 



0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 

(d) Eigenvector 5: A| = 96 





0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 

(e) Eigenvector 7: A^ = 129.5 



0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 

(f) Eigenvector 12: Af 2 = 200.3 



Figure 21: Eigenvectors and corresponding eigenvalues for the domain obtained 
for p xy = 20%, p xz = —10%, p yz = —60%. The mesh is constructed using 1600 
points and is shown in figure 18 
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The eigenfunction expansion for Green's function can be written using the 
previously computed eigenvectors and eigenvalues: 



oo 

G(r, r', = ]T C n g n (r, r')\M^', 6') 



71=1 



The coefficients C„ can be computed by imposing the initial condition for 
Green's function: 

G(0, r', ^, 0') = -r^-H-Stf - r )%' - ^ )5{9' - 6 ). 
Tq sin 6»o 

Since we have ensured the initial condition <? n (0,r') — — r ) for function 

g, we obtain the following equation for the coefficients C n : 



oo 

V C„*„(^, 6') = -t-^SW - <po)6(& - 9 ) 
sin On 



(41) 



em are or- 



n=l 

Given our weak formulation (40 1, the eigenvectors for our probl 
thogonal for the scalar product weighted by sin0': 

/ / 6')* m (<p', 0') sind'dcp'de' = s n , m . 

J Jn 

We multiply equation (41 ) by ^ m ((p', 8') sin 9' and we integrate over the whole 
domain: 

C m = [[ — ^-* ro ft/, 9 1 ) sin 0'%' - - OoWM', 

J Jo. sm6> 

and hence C m = ^ m ((po, Oq). The final formula for Green's function is: 



II , 2 
_ +r OO 

e 2t 



T 



(42) 



6.3 Joint survival probability 

Similarly to the two dimensional case, we denote by Q(t,T,x,y,z) the joint 
survival probability of issuers x, y and 2 to a fixed maturity T. This solves the 
following pricing equation 

Qt ~t~ ^:Qxx ~t~ ~zQyy 7\Q ZZ PxyQxy ~t~ PxzQxz ~t~ PyzQyz 0; 

(43) 
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with final condition Q(T,T, x,y, z) = 1 and zero boundary conditions. We 
proceed to a similar change of variables as described in section 6.1 and using 
the expression for Green's function given in equation (42) we obtain (t = T — t): 

i>oo p-nj p&(tp) 



Q(i~, r Q ,ip ,9 ) = 



G (t, r ,r',ip , <f', do, 0') r' 2 sin 9' dO'dip'dr' 



o Jo Jo 

oo 

n=l 
oo 



0') sin 6' dtp 1 dB 1 



OO r' 2 +rg 

"e 2 T 



i r 2 rfr 



g^-! r(^ + f) 

r(o„+i) lCl 



'«- L y 4. 1 _ UL 



(44) 



x*„((^o,0o) J J^ n (v',O')sin0'dip'de', 

where ±F\ denotes the confluent hypergeometric function and v n - 
We observe that this is a generalization of equation ( 26 ) , which we obtained in 
the two dimensional case. 



6.4 Application to the CVA computation 

We associate the process xt with the protection seller, the process yt with the 
reference name and z t with the protection buyer. The pricing equation for 
computing CVA or DVA in the case where all three names are risky is given by: 

V t + \V XX + ^Vyy + ^V zz + P X yV xy + P XZ V XZ + PygVy z ~ QV = 0, (45) 

with the final condition V(T, T, x, y, z) — and boundary conditions depending 
on the payoff. 

In the case of the CVA calculation, a payout is due if the protection seller 
defaults. If we denote by Rps the recovery of the protection seller, the payout 
is: 

V CVA (t, T, 0, y, z) = (l- R PS ) V(t, T, y)+, (46) 

where V(t, T, y) + is the positive value of the single name default swap with 
non-risky counterparts at the time of the default of the protection seller. 
Similarly we have the payout for the DVA calculation: 

y DVA (;,7>,j/,0) = (1 - R PB )V(t,T,y)-, (47) 

where Rpb is the recovery of the protection buyer and V(t, T, y)~ is the negative 
value of the single name default swap with non-risky counterparts at the time 
of the default of the protection buyer. 

For both CVA and DVA calculations the boundary conditions are for all 
other cases. 



Following the same procedure as in section [67f| the function and first variable 
changes (see equation ( 32 ) ) are applied such that the pricing equation becomes: 

1„ 1„ 1. 



Ut + ^U aa + ^Ufsp + -Ujj 



0. 



with the final condition U(T, T, a, ft, 7) = and boundary conditions except 
for: 

U OVA (t, T, 0, ft, 7) = e e(T - 4) (1 - Rps) V(t, T, p xy ft) + , (48) 
in the case of the CVA calculation, and 



U DVA (t,T,a,ft, 



-PxzP xy a+/3 



} = e« T - t )(l-R PB )V(t,T,p xy a + p aiv 0) , (49) 



for the DVA calculation. 

The second change of variable is applied (see equation ( 33 )) and the modified 
pricing problem is: 



Ut 



1 



1 0* 

r dr' 



(rU) 



1 



1 



1 d 



0, 



(50) 



v sm~ t> sin 9 89 

with final condition U (T, T, r, <p,9) = and boundary conditions except for: 

U CVA (t, T, r, 0, 9) = e ^ T ~*) (1 - R PS ) V(t, T, p xy r sin 9) + , (51) 
for the CVA calculation, and 



U DVA (t, T, r, v , Q(<p))= e^ T -*)(l - R PB )V(t, T, (p xy sin tp + p xy cos <p) r sin 9) . 

(52) 

for the DVA calculation. 

We denote by C the Laplace operator in spherical coordinates: 



CU = 



1 



1 d 2 



1 



2 /) V<ptp 



J_d_ 
\^9d~9 



{smOUe) 



In order to obtain our solution U that satisfies the pricing equation (50) we 
start from the following identity: 



OO /' zu 



t JO JO JO 



e(v>) 



(U t + CU) G (t' - t, r, if, 9) r 2 sin 9 d9 dip dr dt' = 0, (53) 



and perform a scries of integration by parts. As in the two dimensional case, 
we use the boundary conditions, the initial condition for Green's function and 
final condition for U, along with the fact that G t — CG = 0, and we obtain the 
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final pricing formula for U: 



U (t,T, r ,(p ,e ) = 

Tooro 

' sme(ip) U (t 1 , T, r, (p, 9 (<p)) G g (*'- t, r, <p, 9 (</?)) dtpdrdt' 



t o o 

Too oo 



U(t',T, r, tp (w), 6 (w)) G„(f - t, r, y>(w),e (w)) 



t 
cT /■oo /.efTxj 



2 Jt Jo 



sin 9 (w) 

[/ (tf, T, r, w, 9) G v (*' - t, r, w, I 
sin 9 



Q u (u)) diudrdt' 



-dOdrdt' 



1 f T f°°f e ^U(t',T,r,0,0) G„(f -t.r.O, 



v_ 

sin 



-dOdrdt'. 



(54) 



2 Jt Jo 

We note that for one of the integrals above we have used the parametric 
representation of the boundary of our domain given by formulas (33) and (34). 
To obtain the precise formulas for the CVA and DVA calculations we use the 
boundary conditions in equations (51) and (52) respectively: 



Tooe(O) 



[/ CVA (i,7> ,ip ,e ): 



U CVA (t',T,r,0,9)G v (f-t,r,0 , 9) 
sm6 



t o o 



(55) 



U OVA (t,T,r ,<po,do) = 

Too w 

-\JJ JsmO(ip) U DVA (t', T, r, <p, 9 (<p)) G e (t'- t, r, tp, 9 {<p))d V drdt' 
t o o 

+ 1 JJJU^ T, r, m 9(c)) Git'- t, r, 9 H ) 
■> 1 1 1 sin 9 (w) 



t o o 



(56) 



These original formulas provide a new way of consistently computing the CVA 
and DVA. Similar ideas can be used for many other purposes, which will be 
discussed elsewhere. 



7 Numerical results 

In this section we present the results of the CVA and DVA calculations for a 
single name credit default swap. We compare the breakeven coupon obtained 
for a standard CDS to the ones obtained when either the protection buyer or 
the protection seller are risky (using the 2D formulation and results) , as well as 
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when both are risky (using the 3D formulation and results) . When using the 2D 
formulation and considering that either the protection seller or the protection 
buyer are risky, the two parties will not agree on the breakeven coupon of the 
CDS. This problem goes away when using the full three dimensional framework, 
where both are risky, and the problem becomes symmetrical. 

We consider three issuers for our example: X as a protection seller, Y as the 
reference name of the CDS and Z the protection buyer j^J We have chosen risky 
entities for the protection seller and the protection buyer such that the effect of 
the CVA and DVA adjustments on the break even coupon are non negligible. 
We calibrate our inputs to the model to market data from the 15th of December 
2011 (see table [TJ. 



Inputs 


X 


Y 


Z 


Initial value 

(T 

Recovery 


0.0359 
2.44% 
50% 


0.3035 
10.45% 
40% 


0.1199 

6.3% 
40% 



Table 1: Input parameters calibrated to market data (15th December 2011). 



The initial value is a measure of the relative distance to default. This has 
been obtained using the share price on that date, together with the outstanding 



number of shares and total liabilities for that company (see Lipton and Sepp 



2009) for a detailed description of the calibration). The volatility a has been 



calibrated such that the 5Y single name CDS spread is matched to the market 
spread (the 5Y point has been chosen as it is usually the most liquidly traded 
contract). 

For the two and three dimensional cases we also need the correlations be- 
tween the different issuers as inputs to our model. These can be calibrated 
from the prices of first to default swap contracts if such contracts including the 
relevant names are available on the market. Alternatively, we can proxy these 
correlations by assigning a sector to each issuer and then using the sector-to- 
sector historically estimated correlations^] In this section however, we aim to 
show the impact of CVA and DVA on the breakeven spread of a CDS, and 
hence we use different sets of pairwise correlations for the same group of issuers 
in order to illustrate a variety of cases. 

Figure [23] shows the simple case where all the pairwise correlations are 0. 
In this simplified case we can compare our joint survival probability obtained 
through the 3D formulation with simply the product of the individual survival 



probabilities (see figure 22 1. The agreement is very good. 



Figure ^3^ shows the CDS breakeven coupon in different cases. We observe 



3 The issuers chosen for the numerical example are real traded entities and the inputs are 
calibrated to the real market data. 

4 In regulatory capital charge models one needs to estimate sector-to-sector and region- 
to-rcgion correlations. This can be done for example by constructing proxy-portfolios for 
each sector using all the issuers that belong to it and averaging their CDS spreads and then 
computing the correlations of the increments of the time series obtained for different sectors. 
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Joint survival probability 




1.00 2.00 3.00 4.00 5.00 

Time to maturity (years) 



Figure 22: Joint survival probability of the three issuers for p xy = 0%, p xz = 0%, 

Pyz = 0%. 



that the spreads are hyper-exponentially flat at 0, which is a known problem 
of models without jumps. However, for longer maturities we can match well 
against market prices as well as analyse the effect of considering the protection 
seller or the protection buyer or both as being risky. If the protection seller is 
risky, the probability of it non paying the full amount due in the case of the 
default of the reference name is non zero, and hence the protection buyer pays a 
lower coupon as it takes on that risk as well. If the protection buyer is risky, the 
breakeven coupon moves in the opposite direction and the two counterparties 
no longer agree on the coupon. The three dimensional case, where both are 
considered risky, solves this problem as it becomes symmetrical. 
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(b) Change in BEC (compared to the standard CDS 
with non-risky counterparts) 



Figure 23: Impact of counterparty adjustments on the break-even coupon of a 
CDS: p xy = 0%, Pxz = 0%, p yz = 0%. 
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Figure 24 shows the case where the protection seller is highly correlated to 
the reference name. In the case of a default of the reference name, the protection 
seller is likely to default as well, and hence the shortfall between the contractual 
payout and what will actually get paid can be significant. The break-even 
coupon will get adjusted accordingly and will be lower than on a standard fully 
collateralised CDS as the expectation of the payout is lower from the protection 
buyer's point of view. 
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Figure 24: Impact of counterparty adjustments on the break-even coupon of a 
CDS: p xy = 80%, Pxz = 50%, Pyz = 30%. 
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Figure 25 shows the case where the protection buyer is highly correlated 
to the reference name. Since on the default of the reference name the coupon 
payments stop regardless of what happens to the protection buyer, the impact 
of considering the protection buyer as risky in this case is not significant. 
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Figure 25: Impact of counterparty adjustments on the break-even coupon of a 
CDS: Pxy = 20%, Pxz = 30%, Pyz = 80%. 
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Figure 26 shows the case where the protection buyer is highly anti-correlated 
to the reference name. This is intuitively the case where the DVA is largest as 
it is in the cases where the reference name does not default that the protection 
buyer is more likely to default on its coupon paying obligation. This leaves the 
protection seller with a potential shortfall. 
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Figure 26: Impact of counterparty adjustments on the break-even coupon of a 
CDS: Pxy = 20%, Pxz = -10%, Pyz = -60%. 
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8 Conclusion 

This paper contains several useful and original results. First, a 3D extension of 
the structural default framework, where the joint dynamics of the firms' values 
are driven by correlated Brownian motions is proposed. Second, the need for 
such an extension for consistent computation of the CVA and DVA is explained. 
Third, a novel method for obtaining a semi-analytical expression for the Green's 
function combining the eigenvalue expansion technique with the finite clement 
method is developed. As might be expected, in the 3D case, a fully analyti- 
cal expression based on the eigenfunction expansion is not available, since the 
eigenvalues and eigenvectors have to be computed numerically via the finite el- 
ement method. However, given a triplet of correlations, these quantities can 
be precomputed, which allows efficient computations across a range of initial 
points, volatilities or other trade-related data (coupons, recoveries etc.), with- 
out repeating the numerically expensive part. Fourth, it is shown how to use 
the Green's function in order to compute joint survival probabilities for three 
different companies and to calculate the CVA and DVA for a standard CDS. 
Fifth, concrete examples calculating the CVA and DVA for a typical CDS with 
real market data are considered and it is demonstrated that, not surprisingly, 
these adjustments can be very large. It is also shown that only simultaneous 
and consistent computation of the CVA and DVA can explain market clearing 
price for the reference CDS. 
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